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INTRODUCTION 


In  this  report  the  pixel  responses  from  an  indium  antiomide  (InSb)  staring  array  are 
analyzed  for  their  1/f  characteristics.  The  InSb  array  has  a  spectral  band  of  3  to  5.5 
micrometers  (Jim).  Examples  of  the  test  data,  the  tabulated  values,  and  the  camera 
parameters  are  presented  later  in  the  Amber  Camera  Test  Results  section. 

The  basic  hypothesis  presented  in  this  report  is  that  fractional  Brownian  motion  (fBm) 
is  the  proper  model  for  spatial  noise.  The  emphasis  in  this  report  is  on  fBm,  wavelets,  and 
statistics.  The  interested  reader  can  consult  the  reports  of  Scribner  and  others  (References  1 
through  4)  and  the  references  therein  to  find  further  discussion  of  detector  physics  and 
supporting  pcriodograia  spectral  plots  that  clearly  exhibit  a  1/f-typc  characteristic  on  the 
arrays  studied.  Although  the  definitive  work  on  the  physics  of  1/f  noise  has  apparently  not 
been  written,  a  body  of  literature  exists  that  offers  physical  theory  for  the  presence  of  1/f 
noise  in  detectors  (Reference  5). 

The  noise  on  infrared  arrays  has  both  a  temporal  and  a  spatial  component  (Reference 
6).  Spatial  noise  (pixel-to-pixel  variation)  manifests  itself  across  the  array  as  a  nonuniform 
response  by  the  individual  detectors  to  a  uniform  illumination.  Temporal  noise  is  observed 
by  sampling  the  time  series  response  from  a  single  pixel.  According  to  Silverman  and 
others,  temporal  noise  results  from  the  “randomness  inherent  in  the  arrival  and  detection  of 
infrared  photons.”  “The  magnitude  of  temporal  noise  increases  as  the  square  root  of  the 
measurement  time,  whereas  the  magnitude  of  the  signal  increases  linearly”  (Reference  6). 

Scribner  and  others  in  a  series  of  papers  have  investigated  the  1/f  character  of  fixed 
pattern  noise  in  staring  arrays  (References  1  through  4).  According  to  these  investigators, 
“Two  sources  of  spatial  noise  after  nonuniformity  correction  are  1/f  noise  and  pixel 
nonlinearities”  (Reference  2).  For  example,  the  report  exhibits  spectral  plots  for  mercury 
cadmium  telluride  (HgCdTe)  starring  arrays  that  show  a  spectrum  with  a  1/f-type  noise 
characteristic  in  the  lower  frequency  band  and  with  a  white-noise  spectrum  in  the  remaining 
frequency  band.  The  plots  also  show  that  platinum  silicide  (PtSi)  arrays  have  a  spectrum 
that  is  mostly  white.  The  HgCdTe  had  a  frame  rate  of  60  Hertz  with  a  spectral  band  of  4.33 
to  4.9  jim,  and  the  rate  for  the  PtSi  3.6-  to  4.3-  fi m  array  was  30  Hertz.  This  report  is  not 
the  place  for  a  detailed  review;  the  reader  interested  in  the  details  is  referred  to  References  1 
through  4. 

Noise  having  a  power  spectrum  that  is  inversely  proportional  to  the  inverse  of  the 
frequency  (i.e.,/*1)  is  classic  1/f-type  noise  and  is  sometimes  called  pink  noise.  It  has 
equal  power  in  octave  frequency  bands.  Flicker  noise  or  1/f-type  noise  is  often  identified 
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by  a  -1  decade  per  decade  asymptotic  slope,  which  it  produces  on  a  log-log  plot  of  the 
power  spectrum — specifically,  the  periodogram.  Pink  noise  is  encountered  in  a  wide 
variety  of  physical  systems,  including  semiconductor  devices.  The  papers  by  West  and 
Shlesinger  (Reference  7)  and  Montroll  and  Shlesinger  (Reference  8)  present  highly 
readable  introductions  to  the  classic  1/f  noise.  Noise  having  a  power  spectrum  that  is 
inversely  proportional  to  the  inverse  square  of  the  frequency  (i.e.,/'2)  is  called  brown 
noise  in  an  allusion  to  Brownian  motion.  Black  noise  has  a  power  spectrum  that  is 
inversely  proportional  to  some  power — for  example,  a  of  the  frequency — where  the 
power  law  exponent  is  greater  than  two.  Black  noise  is  beyond  brown  noise.  All  these 
different  colored  noises  whose  power  spectrums  are  inversely  proportional  to  some  inverse 
power  a  of  the  frequency  (i.e.,/*a  for  1  <  a  <  3)  are  examples  of  1/f-type  noise 
processes.  White  noise  has  a  flat  spectrum  over  some  finite  range  of  frequencies.  For  a 
lively  and  intriguing  discussion  of  these  colored  noises,  see  the  book  Fractals,  Chaos, 
Power  Laws  (Reference  9). 

Implicit  in  the  discussion  of  1/f  processes  is  the  fact  that  they  have  a  well-defined 
spectrum.  The  theory  of  flicker  noise  is  delicate  and  far  from  trivial,  because  the  frequency 
spectrum/*1  is  not  integrable  on  the  semi-infinite  interval  (0,»),  which  suggests  infinite 
variance  and  therefore  that  the  spectrum  of  the  process  apparently  cannot  exist.  Normally 
these  characterizations  come  from  stochastic  process  theory.  Stochastic  process  theory 
provides  a  dynamic  mathematical  theory  of  random  events  that  occur  in  time  or  space. 
Recently  Solo  has  addressed  the  nonintegrability  issue  of  the  1/f  spectrum;  namely,  he  has 
rigorously  established  “An  explicit  definition  of  the  spectrum  . . .  that  is  compatible  with 
accepted  definitions,  yet  allows  nonintegTable  spectra”  (Reference  10).  His  succinct 
statement  of  the  problem  clearly  implies  that  standard  spectral  theories,  however  useful,  do 
not  address  the  nonintegrable  case.  Fortunately,  his  new  results  include  fBm  with  a 
spectrum  obeying  f-°  for  1  <  a  <>2,  which  is  the  interval  of  interest  for  staring  focal 
plane  arrays  (FPAs).  Fractional  Brownian  motion  noise  can  model  all  of  these  different 
colored  noises.  The  critical  step  in  determining  a  model  for  the  array  noise  is  to  estimate  the 
parameter  a,  something  that  is  hard  to  obtain  from  measured  data.  However,  we  do  have  a 
set  of  candidate  statistics  that  we  will  use  to  estimate  a.  It  should  be  noted  that  Scribner 
and  others  do  not  postulate  a  fractional  Brownian  motion  model  (Reference  1).  However, 
they  do  report  a  mean  pixel  value  fora  of  1.2  for  a  midwave  infrared  FPA. 

Published  reports  clearly  recapitulate  the  existence  of  some  type  of  1/f  noise  in  FPAs, 
and  recent  important  theoretical  progress  clarifies  the  link  between  fractional  Brownian 
motion,  wavelets,  and  periodograms.  If  a  strong  case  can  be  made  for  fBm  as  a  reasonable 
model  for  staring  array  spatial  noise,  then  a  paradigm  for  point  target  detection  in  fBm  can 
be  formulated.  In  fact,  recent  studies  by  Womell  (References  1 1  and  12)  and  Barton  and 
Poor  (Reference  13)  address  the  problem  of  signal  detection  in  fractional  Brownian  noise. 
Their  approach  is  similar  to  that  already  used  with  respect  to  radar,  in  which  a  target- 
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detection-strategy  paradigm  is  defined  by  a  maximum  likelihood  algorithm.  The  algorithm 
defines  an  optimum  target-detection  strategy  with  a  false-alarm  rate  that  must  discriminate 
between  fluctuating  targets  embedded  in  white  noise  and  random  clutter.  Although  such 
methods  are  not  perfect,  they  are  useful  for  design  and  do  define  detection  algorithms. 
Thus,  a  determination  of  the  array  noise  is  an  important  basic  step  in  establishing  a  similar 
detection  agenda  in  staring  FPAs.  Another  outgrowth  of  having  a  model  for  array  noise  is 
that  novel  array  nonuniformity-compensarion  algorithms  that  incorporate  fBm  might  be 
developed.  Nonuniformity-compensation  algorithms  are  necessary  to  correct  for  pixel 
nonuniformities,  but  1/f  noise  and  system  instabilities  foil  many  algorithms;  thus, 
continuing  recalibration  is  required. 

In  the  next  section  of  this  report,  FBm,  Wavelets,  and  ARIMA  Models,  some  important 
properties  of  fBm  and  wavelet  transforms  are  reviewed.  A  key  parameter  in  fBm  is  then 
identified  with  the  exponent  a  in  the  1/f-type  spectra.  Some  concepts  and  operations  from 
stochastic  theory  that  are  relevant  to  the  definition  of  fBm  are  defined.  Next,  some 
statistical  algorithms  are  introduced  that  can  be  used  to  directly  measure  the  parameters  in 
fBm.  Autoregressive  moving  average  models  with  fractional  differencing  first  introduced 
by  Hosking  (Reference  14)  are  reviewed.  Hosking’s  model  provides  simple  and 
fundamental  insight  about  the  kind  of  discrete  stochastic  processes  that  obey  f'a  behavior 
in  a  neighborhood  of  the  origin.  This  discussion  explains  why  fractional  differencing  and 
fractional  sums  are  fundamental  operations  that  are  necessary  to  properly  model  fBm. 

In  the  FBm  Simulation  and  Examples  section  of  this  report,  concepts  and  statistics 
introduced  in  the  FBm,  Wavelets,  and  ARIMA  section  are  illustrated  by  using  a  Hosking 
algorithm  to  synthesize  fBm.  After  the  algorithm  is  introduced,  sample  paths  and 
increments  of  fBm  are  presented  along  with  their  periodograms,  autocorrelations,  wavelet 
transforms,  and  marginal  distributions.  Some  small  sample  studies  of  the  statistics  that 
have  been  proposed  to  estimate  the  fBm  parameters  are  also  presented. 

Some  details  about  the  Amber  camera  experimental  test  conditions,  some  processed 
examples,  and  some  measured  parameters  from  the  test  are  presented  in  the  Amber  Camera 
Test  Results  section.  Some  time  samples  from  selected  individual  pixels  and  their 
increments  are  presented  in  a  multipanel  format  along  with  their  periodograms, 
autocorrelations,  wavelet  transforms,  and  marginal  distributions.  Thus,  the  analysis  in  this 
section  is  exclusively  concerned  with  spatial  noise.  After  the  multipanel  graphs,  a  table  of 
key  fBm  parameters  is  presented  that  summarizes  the  results  of  60  selected  pixels  that  are 
sampled  across  the  array. 
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FBm,  WAVELETS,  AND  ARIMA  MODELS 


A  flicker  noise  (FN)  or  1/f-type  process  is  a  rime  series  with  a  low-frequency  spectrum 
F  (co)  of  the  form  co'a  a  >  0  in  a  neighborhood  of  the  origin.  When  a  =  1  then  the  rime 
series  exhibits  a  long-term  temporal  dependency.  The  study  of  such  processes  is  somewhat 
paradoxical,  because  for  a  =  1  the  spectrum  is  not  integrable  on  the  semi-infinite  interval 
(0,<»).  This  suggests  that  the  model  implies  infinite  variance  and  therefore  that  the 
spectrum  of  the  process  apparently  cannot  exist.  Recently,  Solo  has  shown  that  fractional 
Brownian  motion  is  a  mathematical  model  for  FN  that  has  a  mathematically  correct 
nonintegrable  spectrum  of  the  form  ara  with  1  <  a  <,  2  (Reference  10). 

Fractional  Brownian  motion  is  an  example  of  a  stochastic  process.  Stochastic  process 
theory  belongs  to  the  “dynamic”  part  of  probability  theory  and  provides  a  mathematical 
theory  of  random  events  that  occur  in  time  or  space.  Stochastic  processes  include  models 
for  the  random  behavior  of  stock-market  indexes,  for  the  path  of  a  particle  in  Brownian 
motion,  for  the  number  of  particles  emitted  by  a  radioactive  source,  for  ocean  clutter  fields, 
and  for  the  behavior  of  queues,  to  name  a  few  common  examples.  In  this  report  the 
fractional  Brownian  model  is  proposed  as  a  stochastic  process  that  models  the  detector 
noise  field  present  on  infrared  staring  FPAs.  Moreover,  if  one  subscribes  to  the  theory  that 
1/f  noise  is  the  ubiquitous  spectrum  for  detector  and  semiconductor  noise,  then  fBm  may 
be  expected  to  be  a  relevant  model  for  many  other  devices. 

Because  the  mathematical  details  are  complex  and  are  available  to  the  interested  reader 
in  References  15  and  16,  only  the  highlights  of  the  necessary  stochastic  concepts  are 
presented  here.  In  general  a  stochastic  process  is  a  function  X  (■  ,• )  of  at  least  two 
variables.  One  variable  describes  the  observation  point  t,  and  the  other  variable  represents 
the  random  event  space  O.  The  space  of  random  variables  Q  describes  a  set  of  random 
events  that  coincide  with  the  observation  point.  Generally  the  dependency  of  the  stochastic 
process  X  (•  ,* )  on  the  space  of  random  variables  £2  is  suppressed,  and  only  the 
dependency  on  the  time  variable  is  explicit.  The  observed  time  series  X(t)  represents  a 
“realization”  or  “sample  path”  of  the  stochastic  process  X  (• ,-  );  thus,  the  value  X(t)  at 
each  time  value  is  a  random  variable  that  depends  on  Q .  Generally  the  dependency  on  Qis 
suppressed. 


The  average  or  expected  value  operation  is  of  course  useful  in  describing  the  properties 
of  the  stochastic  process  X(t).  Some  of  the  common  operations  are  denoted  as  follows: 
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mean  E(X(t)),  autocovariance  Cxx  (f)=E(X*(r)X(r))  -  ECX(t))E(X(t)),  where  * 
denotes  complex  conjugate  and  autocorrelation  function  C XX (*)/ yjCxxU)  C  XX  (* ) 

Given  a  set  of  discrete  samples  of  a  stochastic  process,  the  autocorrelation  becomes  the 
sum  for  M  <  N  (Reference  17): 

^  A7— jm( — 1 

cxx{m)~~  lLx{n)x{n+m),  0<m£M  -1 


If  the  sample  mean  is  first  subtracted  from  x(n),  then  is  an  estimate  of  the 

autocovariance,  and  the  autocovariance  divided  by  the  square  root  of  c^Qj  (standard 
deviation)  is  the  correlation.  Sometimes  the  term  “auto*’  is  omitted  when  the  operational 
variables  are  unambiguous. 

The  concept  of  stationary  for  a  time  series  is  analogous  to  the  concept  of  steady-state 
analysis  for  time-invariant  systems.  Roughly  speaking,  a  stochastic  process  is  called 
stationary  if  some  major  property  is  invariant  with  respect  to  a  time  translation  of  the  time 
axis.  A  stochastic  process  defined  on  some  space  Cl  and  over  some  time  interval  T  is 
stationary  if  the  statistics  are  not  affected  by  a  shift  in  the  time  origin.  This  means  that  the 
statistics  of  X  (t)  and  X  ( t  +  x  )  are  the  same  for  all  x  0. 

A  stochastic  process  defined  on  some  space  Cl  and  over  some  time  interval  T  is  weakly 
stationary  or  covariance  stationary  if  the  mean  and  variance  are  the  same  for  all  time  and  if 
its  covariance  cov(X(t),  X(s))  depends  only  on  the  absolute  difference  |  s  - 1 1  for  every  s 
and  t  in  T. 

If  a  process  is  stationary  then  it  is  weakly  stationary;  however,  the  converse  is  not  true, 
because  covariance  stationary  involves  only  the  first  two  moments.  A  stochastic  process 
X(t)  defined  on  some  space  Cl  and  over  some  time  interval  T  has  stationary  increments  if 
the  probability  properties  of  the  process 

VX(r)  =  X(r)  -  X(t  -  x) 
depends  only  on  the  lag  variable  x. 

Mandelbrot  and  Van  Ness  (Reference  18)  introduced  the  class  of  stochastic  processes 
called  fractional  Brownian  motion  (fBm),  which  includes  ordinary  Brownian  motion  as  a 
special  case.  The  process  is  denoted  as  Bn(t)  and  is  a  zero-mean  nonstationary  stochastic 
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process  that  is  indexed  by  a  single  scale  parameter  H.  Here  we  assume  that  BH  (0)  =  0  for 
convenience.  The  fBm  covariance  has  the  following  form  for  0  <  H  <  1 : 

2 

cov(BH(t,s»  sE[BH(t)BH(s)]  =  (jsf"  +  |t,2//  -|t-s|2//) 

When  the  index  H  =  1/2,  the  fBm  is  ordinary  Brownian  motion.  With  the  covariance 
equation,  the  following  inequalities  for  the  covariance  equation  can  be  deduced  for  t  >  0 
and  0  <  H  <  1  using  the  simple  chain  of  inequalities  (for  H  <  1/2  a  similar  chain  can  be 
deduced). 

2  22H  <=>  log2  2  >  2Hlog2  2  «  j  >  H 


>  0 


COV(BH  (t)BH  (-t»  =  E[-BH  (t)BH  (-t)]  =  =  0 


<  0 


if  H  > 
if  H  = 


1 

2 

1 

2 


J  H<  - 
2 


When  for  H  =  1/2  the  past  and  future  increments  are  not  correlated,  otherwise  fBm 
exhibits  persistence  or  antipersistence  as  a  function  of  the  index  R  For  H  >  1/2  persistence 
rules,  which  means  past  increasing  or  decreasing  trends  are  on  the  average  replicated  in  the 
future.  Antipersistence  occurs  when  H  <  1/2;  in  this  case  past  increasing  or  decreasing 
trends  are  anticyclical  in  the  future.  Persistence  is  the  presence  of  significant  dependence 
between  observations  that  are  a  long  span  apart,  whereas  a  common  model  for  time  series 
is  that  observations  separated  by  a  long  time  span  are  nearly  independent.  For  this  reason 
fBm  is  useful  for  modeling  the  time  series  of  phenomena  that  exhibit  persistence  over  a 
wide  range  of  time  scales.  These  concepts  and  an  exact  expression  for  the  covariance  are 
found  in  Reference  19. 

It  follows  from  the  covariance  equation  that  the  variance  of  fBm  satisfies 

2 

var(BH  (t ))  =  E[BH  (t)BH  (t)]  =  <r  \x\2H 
and  diverges  with  increasing  time. 
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Although  fBm  is  not  weakly  stationary,  VBj-j  (t)  has  stationary  increments.  Moreover, 

this  increment  process  is  self-similar  for  any  a  >  0  where  equality  means  equality  in 
distribution  (i.e.,  the  statistics  of  fBm  are  invariant  to  time  scaling). 

H 

Bh  (at)  =  a  Bh  (t) 

This  self-similarity  means  that  a  realization  of  such  a  process  is  a  fractal  curve,  and  its 
fractal  dimension  D  is  D  =  2  -  H  (Reference  20). 

A  specialized  definition  of  fBm  Xati)  is  given  by  the  expression  found  in  Ramanthan 
and  Zxeitouni  (Reference  21).  For  the  index  H  0  <  H  <  1  and  the  zero-mean,  stationary' 
white  Gaussian  noise  random  process  w(t)  with  unit  spectral  density,  the  density-  is 

1  °r  1 1  |tf-l/2  i  |H-l/2\  f.  iM-l/2 

x*(,)=r(^5)[Iv'_Tl  HMIM  wU)dT 

where  T  represents  the  Gamma  function.  The  proper  mathematical  definition  of  the 
spectrum  of  fBm  is  a  delicate  matter.  To  justify  fBm  as  a  model  for  1/f-type  processes  or 
FN  means  that  at  low  frequencies  its  spectrum  approximates  co  a.  Recently,  some 
important  new  results  by  Flandrin  (References  22  and  23),  Hosking  (Reference  14),  Solo 
(Reference  10),  and  Womell  (References  1 1  and  12)  have  all  made  significant  progress  in 
refining  the  mathematical  connections  between  fBm  and  its  1/f-type  spectrum.  Hosking 
shows  that  a  family  of  autoregressive  integrated  moving  average  (ARIMA)  processes  with 
fractional  differencing  has  a  co  spectrum.  Flandrin  uses  a  wavelet  transform  to  show  how 
the  spectrum  of  fBm  can  be  estimated  and  defined,  and  Womell  uses  a  wavelet  transform  to 
produce  models  that  approximate  a  1/f-type  spectrum.  Solo  justifies  the  estimation  of  a  in 
fBm  spectra  co  * a  from  n  discrete  time  data  samples  by  the  slope  of  a  log-log  plot.  The  plot 
is  the  log  of  the  periodogram  versus  the  log  of  co.  The  periodogram  is  the  squared 
modulus  of  the  discrete  Fourier  transform  of  the  data. 

To  demonstrate  that  the  spectrum  of  fBm  behaves  as  co  a,  Flandrin  starts  with  the 
wavelet  of  fBm,  which  is  defined  as 

W(b,  a)=-ij  BH  (t)g(— )  dt 
’  *  —80  ^ 


9 


NAWCWPNS  TP  8103 


Here  g  is  any  wavelet  with  scale  parameter  a  whose  Fourier  transform  3  satisfies  the 
admissible  condition 

7  gii-tM2  da,  <  ,  (!) 

-oc  N 

where  3  (g(0)»  is  the  Fourier  transform  of  g(t).  If  3  (g(co))  is  differentiable,  then  it 

oc 

suffices  that  g  be  zero  mean,  i.e.,  jg  (t)  dt  =0  .  for  Equation  2  to  be  satisfied 

—  oo 

(References  24  through  26). 

t  -  $ 

As  shown  by  Flandrin,  the  covariance  of  W(b,a)  is  a  function  only  of - and  thus 

a 

W(t,a)  is  a  weakly  stationary  stochastic  process.  In  other  words,  when  the  nonstationary 
process  BH  (t)  is  analyzed  at  a  fixed  scale,  it  is  weakly  stationary!  Thus,  the  Fourier 
transform  of  cov  (W  (t,a)W(s,a))  is  well  defined.  So  taking  the  average  of  the  Fourier 
transform  3  of  cov  (W  (t,a)W(t,a))  and  integrating  over  all  scales  using  the  normalizing 
measure  yields 


7  3(cov  (W  (t,a)W(t,a))-~y  =  -  (2) 

Thus,  Flandrin  has  established  a  direct  connection  between  the  wavelet  approach  and  the 
1/f-type  process.  Another  remarkable  result  about  the  wavelet  transform  and  fractional 
Brownian  motion  has  been  established  by  Flandrin.  To  explain  his  result  some  additional 
facts  about  wavelet  transforms  are  now  introduced. 

Wavelets  are  not  arbitrary  functions,  and  as  a  minimum  they  satisfy  the  Fourier 
condition  defined  by  Equation  1.  The  wavelets  y/{t)  introduced  by  Daubechies  (Reference 
24)  are  obtained  from  the  scaling  function  <f>(t),  which  is  the  solution  of  the  two  scale 
difference  equation  (a  dilation  equation): 

0(O  =  Ic*4K2f-Jfc) 

i 

If  the  coefficient  sequence  {c*}  is  finite  in  length,  the  wavelet  y/(t)  has  compact  support. 
The  sequence  {c*}  satisfies  the  following  conditions,  which  are  respectively  a  normalizing 
condition,  an  orthogonality  condition  ( $*»  is  the  discrete  Dirac  delta  function),  and  a 
regularity  or  vanishing  moment  condition. 
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Ic*  =  V2,  Ickc*.2w  =  25o«,  I(-Dkr*o  =  0,  m*0,l,  -.p-l 

The  wavelet  equation  is  constructed  from  the  sequence  {c*}  by  the  quadrature  mirror 
condition 

V'CO  =  £(-!)*  c,-*0(2r-l) 


The  Haar  wavelet  has  compact  support  is  quite  simple  and  is  defined  as  follows: 


f+1 


Vfir)  = 


-1 

0 


0  <  t  <  1/2 
1/2  <  t  <  1 
otherwise 


The  scaling  function  is 


4>(t)  = 


0  <  t  <  1 
otherwise 


Although  Haar  wavelets  are  simple  and  compact  in  the  time  domain,  they  are  not  well 
localized  in  Fourier  space  where  they  decay  very  slowly  as  a)'1.  Wavelets  that  are  well 
localized  in  both  time  and  Fourier  space  have  been  discovered  only  recently  by  Daubechies 
(Reference  24).  Wavelets  that  decay  as  rapidly  as  (o  in  Fourier  space  and  exponentially  in 
the  time  domain  are  discussed  by  Mallat  (Reference  27).  Wavelets  that  are  well  localized  in 
both  time  and  space  provide  a  reasonably  sharp  decomposition  of  the  signal  in  both 
domains.  The  localization  is  limited  to  bands  of  nonzero  thickness  in  both  t  and  co  space 
because  of  the  Fourier  or  Heisenberg  uncertainty  principle. 

The  Hilbert  space  of  measurable  square  integrable  functions  is  denoted  by  L2(R) ,  i.e., 
the  space  of  signals  fit)  with  finite  energy,  satisfying 


7|/(0|2d'<~ 


The  compactly  supported  Daubechies  wavelet  forms  an  orthonormal  basis  for  L2(R), 
which  means  that  any  function  that  is  a  member  of  L2(R)  can  be  expanded  into  a  wavelet 
series  instead  of  a  Fourier  series.  A  function  that  is  a  member  of  L2(R)  has  “finite  power” 

and  can  be  expanded  into  a  series  using  the  sequence  of  wavelets  |V2^(/(2y/  -  m)J  as  the 
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basis  functions.  Moreover,  for  some  wavelets  all  translates  and  dilates  of  y/(i)  are 
mutually  orthogonal  using  the  inner  product  {*) . 

(V?  y(2kt -  /),  V27 m))  =  5 *  6* 

The  discrete  wavelet  coefficients  that  represent  the  expansion  in  the  scaling  functions 
are  called  the  approximation  coefficients  a  j(n) 

aj(n)  =  -  n)j  =  '&IToof(0<P(2jt-n)dt 

where  f(t)  is  any  function  in  L2(R)  ■  The  wavelet  expansion  coefficients  are  called  the  detail 
coefficients  dj(n),  because  they  represent  the  difference  between  two  successive 

approximations  and  are  computed  as  follows  (Reference  25). 

di  (*)  =  (f(0, V2V(2'f-  /»))  =  &IZof(0Y(2jt-n)di 

For  any  finite  integer  J,  {-^2*  y/(2J  t — m  )UV?V(2'f  -m)J  is  an  orthonormal  basis  for 
L2(R).  Thus,  any  function  f(t)  in  L2(R)  can  be  written  as  the  infinite  series 

/(r)=V 2J  Xa(m)0(2^t-m)+  X  I  ^2J d j(m)yK2J t ~ m) 

m=—°°  j=Jm=—°° 


This  series  expansion  is  called  a  multiresolution  decomposition  of f(t).  The  decomposition 
is  multiresolution,  because  the  approximation  coefficients  and  the  detail  coefficients  are 
computed  at  successively  different  scales  indexed  by  j.  Note  the  different  role  played  by  the 
approximation  coefficients  and  the  detail  coefficients  in  the  expansion.  Of  course,  for  any 
discrete  sequence  {*<>»•  the  discrete  wavelet  transform  is  actually  computed.  A 

detailed  and  insightful  discussion  of  the  discrete  wavelet  transform  is  found  in  the  report  by 
Shensa  (Reference  28). 

The  diagram  in  Figure  1  illustrates  the  standard  orthonormal  fast  wavelet  transform. 
The  solid  vertical  line  represents  the  division  between  the  analysis  (decomposition)  and 
synthesis  (reconstruction)  parts  of  the  transform.  The  input  vector  x  represents  a  discrete 
data  sequence  that  is  convolved  on  the  parallel  branches  via  the  two  discrete  filters  g  and  h. 
The  coefficients  for  the  two  filters  are  obtained  from  the  wavelet  and  scaling  function. 
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respectively.  The  arrow  i  2  means  drop  every  second  sample.  The  2  is  the  dilation  factor 
for  use  with  wavelets  with  one  octave  bandwidth.  The  output  from  the  top  branch 
represents  the  detail  coefficients  dj(n),  and  the  output  from  the  bottom  branch  represents  the 
approximation  coefficients  aj(n).  If  the  wavelet  transform  uses  orthogonal  wavelets,  then 

either  or  both  sets  of  output  coefficients  may  be  used  as  the  new  input  vector  x.  When  both 
dj(n)  and  aj(n)  are  used  as  separate  and  parallel  inputs  for  the  analysis  stage,  then  a  wavelet 

packet  hierarchical  or  tree-structured  signal  decomposition  is  produced.  The  synthesis  is 

> 

also  hierarchical  in  nature  and  reconstructs  a  final  signal  x  by  using  the  detail  coefficients 
dj(n)  on  the  top  branch  and  the  approximation  coefficients  aj(n)  on  the  bottom  branch.  The 

A  A 

T  2  arrow  means  insert  one  zero  between  samples.  Generally  the  filters  g  and  h  are  related 
to  g  and  h,  but  as  long  as  they  represent  a  valid  orthogonal  wavelet  scaling  function  pair  the 
synthesis  branch  can  use  a  different  wavelet.  The  diagram  in  Figure  1  represents  a  fast 
dyadic  (or  two-band)  wavelet  transform  for  doing  multiresolution  processing  as  introduced 
by  Mallat  (Reference  28)  and  Daubechies  (Reference  24).  Recently,  Zou  and  Tewfik 
(Reference  29)  have  constructed  M-band  generalizations  that  have  a  similar  block  diagram 
with  multiple  parallel  branches  for  the  additional  M-l  wavelet  filters,  and  the  i-2  block  is 
replaced  with  X .  M.  Cohen  and  Daubechies  (Reference  30)  have  shown  how  to  construct 
rational  band  orthogonal  wavelet  bases  with  1 /2-octave  or  even  smaller  bandwidth 
wavelets. 


ANALYSIS 


SYNTHESIS 


> 

x 


FIGURE  I.  Fast  Wavelet  Transform  Block  Diagram. 

Flandrin  has  shown  that  the  detail  coefficients  have  stationary  properties  after  a  wavelet 
transform  is  applied  to  fBm  (Reference  23)!  Tewfik  and  Kim  have  shown  that  the 
vanishing  moment  property  of  the  Daubechies  wavelets  controls  the  correlation  at  a  given 
scale  and  between  scales  (Reference  31).  In  integral  terms  the  wavelet  y/(t)  has  M 
vanishing  moments  if 
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ltm\f/(t)dt  =  0  m  =  0,l,---,M-l  forM  >0 


Not  all  wavelets  have  compact  support,  nor  are  they  all  orthonormal  functions.  In  fact, 
Mallat  and  Zhong  have  shown  that  multiresolution  edge  detection  wavelets  are  not 
orthogonal  (Reference  32).  All  of  these  topics  and  formulas  for  many  other  interesting 
wavelets  can  be  found  in  the  excellent  monograph  by  Daubechies  (Reference  25). 

Hosking  models  1/f-type  processes  by  starting  with  an  ARIMA(0,d,0)  process  with 

1  j . 

\d\  <  —  for  the  0)za  spectrum  and  then  extending  it  to  a  ARIMA(p.d.q)  process 
2 

(Reference  33).  When  d  is  an  integer,  then  ARIMA(p,d,q)  denotes  the  usual  ARIMA 

process  with  p  coefficients  for  the  autoregressive  process,  d**1  integral  differences,  and  q 
coefficients  for  the  moving  average  process.  To  summarize  the  important  properties  of  the 
ARIMA(0,d,0)  process  model  for  fBm,  let  us  define  the  differencing  operator  V,  which, 
when  applied  to  a  discrete  sequence  (discrete  time  signal)  x(n),  yields  the  derived  sequence 

Vx(n)  =  (l-B)x(n)  =  x(n)  -  x(n-l) 


where  B  is  the  backward  shift  operator  Bx(n)  =  x(n-l).  Fractional  differencing  is  obtained 
by  repeated  applications  of  B  d  times,  where  d  is  a  real  number.  The  operation  is  formally 
defined  by  the  binomial  expansion  of  the  operator  V4: 


Vdx{n)  =  (1  -  B)dx(n)  =  ^j(-B)A:x(n) 


=  x(n)  -  dx(n  -1)  -  ^  d(l-d)x(n  -2)  -  \  d(l-d)(2-d)*(n  -3)  - 
2  6 


(3) 


For  \d\  <  ~  the  infinite  series  converges  in  mean  square  and  so  the  operation  V4  is  well 

defined.  The  ARIMA(0,d,0)  is  the  discrete  time  sequence  x(n)  that  is  represented  by  the 
following  equation: 

Vdx(n)  =  e(n) 


where  the  operator  Vd  is  defined  in  Equation  3  and  the  white-noise  Gaussian  input  e(n)  is  a 
sequence  of  independent  and  identically  distributed  random  zero- mean  variables  satisfying 
the  equations 

E(e(n))  =  0 
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E(e(n)e(j))  = 


for  n=j 
for  n*j 


The  Hosking  ARIMA(0,d,0)  model  is  discrete  fractional  Gaussian  noise  (fGn).  Discrete 

fBm  is  then  generated  by  adding  up  the  increments.  The  ARIMA(0,|,0)  process  is  fGn 

with  a  1/f  spectrum,  while  the  ARIMA(0,0,0)  process  is  fGn  with  constant  spectra) 
density.  It  has  some  very  interesting  properties  that  qualify  it  as  a  candidate  for  fBm.  If 

x(n)  is  an  ARIMA(0,d,0)  process  for  \d |  S  -t  then  x(n)  is  stationary  with  1/f-type 

spectral  density: 


S(co)sspectrum((o)  =  (2sm^Q)r2d  0  <co<n  S(a>)-af2d  as  tu->0 

The  correlation  function  of  x(n)  is  as  k  — » «» 

„k_d{\  +d)-(k  -1  +  d)  Td-d),^., 

P  (1  -d)(2  -d)-(k  -d)~  T(d)  k  (4) 

where  T  represents  the  Gamma  function.  The  partial  correlations  of  x(n)  are 

*"=%-  d) 

For  the  d  £  0  ARIMA(0,d,0)  process,  the  correlation  decays  hyperbolically  as  the  lag 
increases,  as  opposed  to  exponential  decay  of  short-term  memory  processes  that  are  typical 
of  autoregressive  moving  average  (ARMA)  (p,q)  processes,  which  are  now  defined. 

In  the  digital  signal  processing  literature,  ARMA  processes  arc  simply  time-invariant 
constant  coefficient  difference  equations.  An  autoregressive  process  (AR(p))  is  a  linear 
difference  equation  with  constant  coefficients  driven  by  white  noise.  The  autoregressive 
process  is  defined  as  a  real  random  sequence  X(n)  with  real  coefficients  iff  _  •  and  input 
e(n): 

X(n)  -  £  VPjX(n-i)  +  e(n)  y/  *  0 

i«l  r'r 

The  AR(1)  model  is  the  first-order  difference  equation 
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X(n)=Yi,iX(n- 1)  +  e(n) 


which  is  stationary  and  has  a  finite  variance  when  jV^i  i|  <  1  •  Moreover,  the  difference 
equation  is  stable  and  has  an  autocorrelation  function  that  decays  exponentially: 

pm=C2cdm)  =  lm 

cxx(0)  ' 

A  comparison  of  the  two  autocorrelation  functions  for  the  AR(1)  and  the 
ARIMA(0,d,0)  processes  clearly  reveals  the  difference  between  persistence  and  short 
memory  processes.  The  autocorrelation  coefficients  for  the  AR(1)  decay  exponentially  and 
are  positive  or  oscillate  in  sign  as  m  — »  «>,  whereas  the  coefficients  in  Equation  4  for  the 
ARIMA(0,d,0)  process  decay  hyperbolically  as  k  -*  °o . 

A  moving  average  process  (MA(q))  is  a  linear  difference  equation  with  constant 
coefficients  driven  by  white  noise.  It  is  defined  as  a  real  random  sequence  X(n)  with  real 
coefficients  tV,  and  white-noise  input  e(n): 


X(n)=l  tfq4e(n-i)+  e(n)  where  *  0 
i=l 

The  ARIMA(0,d,0)  process  has  three  free  parameters:  the  fractional  difference  d,  the 
mean,  and  the  variance  of  the  process.  The  ARIMA(p,0,q)  processes  have  the  p  +  q 
coefficients  as  free  parameters.  Thus,  the  Hosking  ARIMA(p,d,q)  process,  with  p 
coefficients  for  the  autoregressive  process,  fractional  differences  d,  and  q  coefficients  for 
the  moving  average  process,  is  a  combination  of  the  above  models.  The  model  has  the 
ARIMA(p,0,q)  process  to  account  for  the  short-term  properties  and  the  ARIMA(0,d,0) 
process  to  account  for  the  long-term  persistent  behavior  of  the  observed  time  series.  An 
algorithm  for  simulating  both  the  ARIMA(p,d,q)  and  ARIMA(0,d,0)  processes  is  defined 
by  Hosking  (Reference  33)  and  is  used  in  this  report  to  model  these  processes. 

The  ARIMA(p,d,q)  processes  satisfy  the  following  equation  with  mean  fj.  and  fractional 
difference  operator  Vd  of  order  d.  In  this  equation,  e(n)  is  the  white-noise  sequence 
defined  for  the  ARIMA(0,d,0)  process  polynomial  operators  defined  for  the  operator  B . 

<f>(B)Vd(y(n)  -jO=0(B)e<n) 

0(B)sl-<plB-  •••  -<ppBP 
■&{B)  h  1-  $\B-  ••• 
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In  order  to  explain  the  connection  between  the  range  of  the  parameter  H  in  Equation  2  and 
the  d  in  the  fGn  ARIMA(0,d,0)  model,  we  need  to  know  the  relationship  between  the 
spectrum  of  the  image  sequence  x(n)  and  the  spectrum  5,<«)(g>)  of  the  image 

sequence  y(n) 

y(n)=Vdx(n ) 

that  is  obtained  by  applying  the  fractional  differencing  operator  Vd  to  x(n)  d  times.  Based 
on  filtering  considerations  (Reference  34),  the  relationship  between  the  two  spectrums  can 
be  formally  written  as 


where  z=e~‘“.  This  transfer  function  relationship  follows  when  we  note  that  the  z- 
transform  for  (1-V)  is  1-z.  When  the  sequence  is  differenced  once  (i.e.,  d=l),  the 

spectrum  is  multiplied  by  |l-z|2  =  2(1  -cos at).  Next  we  note  that  (1-cos to)  =  2(sin^)2, 
and  we  have  in  a  neighborhood  of  the  origin 

Sy(n)(a)=(-^)S  ,  (<0) 

O)2  V1Jt(n) 


The  equation  H  =  d  +  ^  combined  with  Equation  5  implies  that  the  fGn  ARIMA(0,d,0) 

mm 

model  with  power  spectrum  of  the  form  in  a  neighborhood  of  the  origin  produces  a 

power  spectrum  for  the  fBm  x(n)  of  the  form  V  2H+1.  Thus,  d  in  the  interval  0  <  d  <  — 

/  co  2 

yields  -£//<l,  while  the  interval  ~-<d<  0  yields  0<f/<  — .  This  relationship 
2  2  2 

implies  that  in  the  Brownian  increment  process  with  persistent  behavior,  as  discussed  in 
Reference  19,  0  <  d  <  produces  a  spectrum  for  the  fBm  sample  path  of  the  form  or2”'1 


in  a  neighborhood  of  the  origin  2  <  2H  + 1  <  3,  while  antipersistent  behavior 


-~<  d  <0 
2 


yields  1  <  2H  + 1  <  2.  Thus,  the  Hosking  model  provides  a  general  approach  to  the  discrete 
1/f-type  processes  with  a  simple  and  logical  relation  between  the  parameters.  For 
antipersistent  behavior.  Solo  develops  the  connection  between  FN  and  fBm  by  providing  a 
firm  theoretical  foundation  for  the  proper  interpretation  of  the  power  spectrum  (Reference 
10).  A  paper  by  Haslett  and  Raftery  (Reference  35)  is  a  comprehensive  study  of  time  series 
with  persistent  behavior. 
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One  of  the  fascinating  historical  aspects  of  fBm  is  the  Hurst  rescaled  range  statistic.  In 
fact,  it  is  the  reason  that  Mandelbrot  used  H  for  the  exponent  in  the  power  spectrum  of  the 
sample  path  of  B^U)  (Reference  20).  The  Hurst  coefficient  and  the  d  in  Hosking’s 

ARIMA(0,d,0)  process  are  related  by  the  equation  H=  d  +  1/2.  However,  the  power 

spectrum  — for  the  fGn  ARIMA(0,d,0)  process  applies  to  the  increments 
co 

VB^(t)  =  B^(t  )  -  B^(t- 1).  The  distinction  between  estimates  that  is  applied  to  the 
sample  path  of  a  time  series  (e.g „B^(t))  or  to  the  increments  of  the  time  series 
(e.g.,  VS^(O)  is  important,  and  one  must  be  alert  to  these  alternatives  when  consulting  the 
literature  in  this  field. 

To  compute  the  Hurst  coefficient,  a  statistic  for  the  range  R  of  a  discrete  data  sequence 
is  defined.  An  explicit  expression  for  R  is 

R(t)=  max  X(t,T)  -  min  X(l,T ) 
l<;/£r  lstsr 

where  t  is  a  sampled  value  of  the  time  series  X(t)  and  T  is  the  time  span  considered.  To 
devise  a  statistic  that  could  discriminate  between  long-  and  short-term  time  scales,  Hurst 
uses  the  dimensionless  ratio  R/S  where  S  is  the  estimated  standard  deviation  defined  as 
follows  for  the  time  series  x(t): 


T/= l 


S  = 


^,t  I*®  -  Urf 


J. 

2 


Finally  the  coefficient  H  is  found  using  a  log-log  plot  from  the  empirical  relation: 


Schroeder  reports  that  the  Hurst  coefficient  for  the  Nile  river  level  minima  taken 
between  the  years  622  B.C.  and  1469  B.C.  has  H  =  0.9,  or  a  spectrum  of  a>~2t 
(Reference  9).  The  equation  H  =  d  +  0.5  yields  d  =  0.4  for  the  corresponding  Hosking 
ARIMA(0,d,0 )  process.  The  Hurst  coefficient  thus  quantifies  the  legendary  cycles  of  the 
Nile’s  annual  floods. 
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FBm  SIMULATION  AND  EXAMPLES 


In  this  section  the  ARIMA(0,d,0)  algorithm  of  Hosking  is  used  to  simulate  fBm  with 

spectrum  — -7—  and  fGn  with  spectrum  —7-.  A  number  of  different  statistics  for 
co2H+l  tu2d 

estimating  the  parameters  of  fBm  are  illustrated.  These  statistics  include  the  Hurst 
coefficient,  the  wavelet  transform,  the  periodogram,  and  some  new  iterative  algorithms  by 
Deriche  and  Tewfik  (Reference  36).  The  new  iterative  algorithms  are  used  to  estimate  the 
parameters  in  the  ARIMA(0,d,0)  models.  Throughout  this  section  the  distinction  between 
the  sample  path  B^(t)  of  fBm  and  the  increments  Vfi^(r)will  be  observed. 

An  algorithm  to  simulate  ARIMA(p,d,q)  models  (p  and  q  non-negative  integers, 
1^1  <  tjO  is  derived  by  Hosking  (Reference  14).  When  p  =  q  =  0,  the  algorithm  simulates 

the  sample  increments  of of  size  n  with  a  normal  marginal 

distribution  and  correlation  function  pt.  The  sample  path  B^{t)  is  the  sum  over  VB^(t). 

For  completeness  the  algorithm  for  the  general  ARIMA(0,d,0)  model  is  reproduced  here 
and  defined  as  follows: 

1 .  Generate  the  starting  value  Xo  fr°m  a  normal  or  Gaussian  distribution  N(0,<7o) 
with  mean  zero  and  variance  cr0  ( <7o  =var(  [xa  ,  •  Xm-i} ))  • 

2 .  For  t=l ,  •  •  < ,n- 1 ,  define  the  recursive  equations  for  0tj  j=l,--.,t 

*.=‘ </,-d 

~  ~  K  j  ~  i»"*>t  •  1 

Calculate  m,  =  x,-j  and  v,  =  (1-  <p2a)v,-\.  Generate  the  next  x,  fromA'(m,,v,). 

/•1 

The  multipanel  graphics  in  this  report  were  generated  using  S-PLUS ®!,  which  was 
also  a  valuable  source  for  many  of  the  basic  statistical  algorithms  utilized  in  the  preparation 
of  this  report. 


1  S-PLUS ®  is  a  registered  trademark  of  Statistical  Sciences,  Inc.,  Seattle,  Washington. 
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All  of  the  sample  paths  R^(r)  or  realizations  plotted  in  Figures  2  and  3  for  fBm  start 
with  n  =  3,000  sample  increments  for  VB^(r).  The  sample  paths  B^(t)  or 

realizations  are  generated  by  the  Hosking  ARIMA(0,d,0)  model  (Reference  33).  In  Figure 
2  three  examples  of  simulated  fBm  are  illustrated.  The  sample  paths  are  illustrated  in  panels 
a,  c,  and  e  with  corresponding  d  values  of  0.4,  0,  and  -0.4.  Each  fBm  realization  has  a 
Gaussian  marginal  distribution  with  zero  mean  and  unit  variance.  In  Figure  14  the 
Gaussian  shape  of  the  marginal  distributions  will  be  examined  by  the  use  of  quantile- 
quantile  (Q-Q)  plots  (Reference  37),  and  in  Figure  15  the  shape  of  the  marginal 
distributions  will  be  estimated  by  histograms  and  by  empirical  density  approximation. 
Panels  b,  d,  and  f  are  the  corresponding  periodograms,  which  are  log-log  plots  of  the 
periodogram  X(ru)  versus  the  log  of  the  frequency  variable  6) .  The  periodogram  X(  a>)  is 
the  modulus  of  the  Fourier  transform  3(x(/t))  of  the  discrete  sequence 

X(tu)  =  ^3(x(n))3(x(n))‘ 

where  N  is  the  total  number  of  samples  and  (  )*  is  the  complex  conjugate.  The  index  H  can 
be  estimated  from  the  slope  of  the  dotted  line.  It  is  the  least  squares  fit  to  the  lower 
frequencies  on  the  log-log  plot  of  the  periodogram. 

Tables  1  and  2  compare  the  performance  of  all  of  the  statistics  that  have  been  defined  in 
this  report.  These  include  the  log-log  plots  whose  slopes  are  used  to  estimate  H .  The  one 
new  estimate  for  H  is  by  Raftery  (Reference  35)  and  it  is  obtained  by  using  his  Fortran 
program  called  fracdiff2. 

These  tabulated  results  are  based  on  a  limited  sample  and  are  meant  to  suggest  only  the 
comparative  variability  among  the  different  statistics.  All  of  the  statistics  in  the  first  four 
rows  of  Table  1  arc  based  on  samples  that  used  the  same  value  of  d,  namely,  d  «  0.4  or 
equivalently  H  =  0.9.  The  statistics  in  the  last  row  of  Table  1  are  based  on  samples  that 
used  d  =  -0.4  or  equivalently  H  =  0.1.  The  column  headings  Hurst,  Wavelet,  Spectrum, 
and  Variance  are  defined  by  matching  them  with  Table  2.  All  of  the  former  statistics  are 

A 

based  on  sample  paths  B//,  and  fracdiff,  which  estimates  d  is  based  on  the  increments 
VBh  •  The  Daubechies  wavelet  Dio  was  used  exclusively  in  the  wavelet  estimates.  When 
two  estimates  are  given  in  the  table,  it  is  because  the  statistics  are  based  on  linear  fits 
obtained  by  the  S-PLUS  routines  least  squares  fit  (Ism)  and  least  median  squares 
regression  (1ms). 


2  fracdiff  is  a  function  in  StatLib  obtainable  via  an  anonymous  FTP. 
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TABLE  1.  Comparative  Statistics  for  the  Estimation  of  the  Index  H. 


No.  of 
samples 

fracdiff 
d  (d=0.4) 

Hurst 

H  (H=0.9) 

Wavelet 

H  (H=0.9) 

Spectrum 

H  (H=0.9) 

Variance 

H  (H=0.9) 

1000 

0.412429 

0.81981 

0.850106(ls) 

0.897877 

0.6638(ls) 

0.86927(lm) 

0.7184(lm) 

2000 

0.422846 

0.8529 

0.849(ls) 

0.9078(ls) 

0.7723(ls) 

0.939(lm) 

0.7881(lm) 

3000 

0.40855 

0.82313 

0.8354(ls) 

0.83944 

0.7 1232(ls) 

0.898  l(lm) 

0.79531(lm) 

5000 

0.423665 

0.92279 

0.7964(ls) 

0.84638 

0.8550(ls) 

0.8985(lm) 

0.8883(lm) 

H  (H=0.1) 

H  (H=0.1) 

H  (H=0.1) 

H  (H=0.1) 

5000 

— 

0.2652(ls) 

0.07068(ls) 

0.08 11  (Is) 

0.1347(ls) 

0.2749(lm) 

0.072 12(lm) 

0.0794(lm) 

0.1919(lm) 

As  expected,  the  estimates  for  H  obtained  by  firacdiff  are  the  most  accurate,  because 
the  program  is  based  on  more  comprehensive  statistical  theories.  Unfortunately,  the  version 

A 

used  for  this  study  could  only  be  used  estimate  d  for  the  persistent  case  (i.e.,  0  <  d  <  1/2). 
None  of  the  other  statistics  did  an  outstanding  job.  They  all  give  rough  estimates,  however, 
after  rounding  to  one  significant  figure.  H  is  bounded  in  the  interval  defined  by  the 
maximum  and  minimum  among  the  three  statistics.  Hum,  wavelet,  and  spectrum. 

Table  2  represents  a  small  sample  study  of  the  effects  of  white  noise  on  the  various 
statistics.  These  are  all  estimated  from  the  sample  path  of  an  ARIMA(0,d,0)  model  plus 
additive  white  Gaussian  noise  N(0,4). 

Table  2  clearly  demonstrates  that  the  statistics  are  sensitive  to  the  white  noise;  however, 
this  is  to  be  expected,  because  their  robustness  has  not  been  assessed.  Based  on  this  small 
sample  study,  it  appears  that  Spectrum  is  the  most  reliable  estimator.  This  is  partially  based 
on  the  fact  that  an  inspection  of  the  periodogram  can  be  made  before  a  frequency  band  for 
the  least  squares  fit  is  selected.  The  inspection  procedure  avoids  the  white-noise  knee  in  the 
periodogram  that  is  clearly  evident  by  comparing  Figures  2  and  13.  Clearly,  a  more 
accurate  estimator  like  fracdiff  is  needed.  The  method  of  Womell  (Reference  12),  which 
uses  the  maximum  likelihood  algorithm,  will  undoubtedly  provide  better  estimates  for  the 
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fBm  parameters  H  and  a2  and  for  the  corrupting  white-noise  variance  ow  than  the 
statistics  used  in  this  exploratory  study. 


TABLE  2.  Comparative  Statistics  for  the  Estimation  of  the  Index  H. 


No.  of 
samples 

Hurst 

H  (H=0.9) 

Wavelet 

H  (H=0.9) 

Spectrum 

H  (H=0.9) 

Variance 

H  (H=0.9) 

3000 

0.72244(ls) 

0.65170(lin) 

0.6878534 

0.92753(h) 

0. /7685(lm) 

5925(h) 

0. 63523(lm) 

5000 

0.67 198(ls) 
0.72429(lm) 

H  (H=0.1) 

0.75177(h) 

0.825170m) 

H  (H=0.1) 

0.91378(h) 

0.66926(lm) 

H  (H=0.1) 

0.626050s) 
0.67901  (lm) 

H  (H=0. 1) 

3000 

0.13528(ls) 

0. 194905(lm) 

0.28619(ls) 

0.25856(lm) 

0.03339(h) 

0.03987(lm) 

0.05932(h) 

0.0485780m) 

5000 

0.1745(lm) 

0.1730(ls) 

0.552000s) 

0.7107(lm) 

0.0967(h) 

0.0568(lm) 

0.04815(h) 

0.5684(lm) 

In  Figure  3  the  identical  B^{t)  realizations  from  Figure  2  are  plotted  in  panels  a,  c,  and 

e.  Panels  b,  d,  and  f  in  Figure  3  are  the  corresponding  correlation.  Note  how  panels  b,  d, 
and  f  all  clearly  exhibit  slow'  decay  with  hyperbolic  shape  and  positive  correlation 
coefficients  that  distinguish  persistent  (d  >  0)  and  an tipersi stent  (d  <  0)  long-term  behavior. 
The  mathematical  details  can  be  found  in  Reference  33  and  in  Equation  4.  These 
correlations  indicate  a  significant  dependence  between  observations  well  separated  in  time. 
This  behavior  is  in  contrast  with  observations  that  are  nearly  independent  over  long  time 
spans. 

In  Figure  4  panels  b,  d,  and  f  are  the  corresponding  correlation;,  for  the 
increments  V£^(r)  that  are  plotted  in  panels  a,  c,  and  e.  Even  in  the  increments  the 

hyperbolic  shape  of  the  correlation  function  is  evident  in  panel  b.  As  expected,  a  significant 
change  in  correlation  has  occurred  for  the  fGn  ARIMA(0,d,0)  model  for  d  =  0  and  for  d  = 
-0.4,  and  it  suggests  a  significant  decorrelation. 


Both  Haar’s  and  Daubechies’  scaling  functions  and  wavelets  are  plotted  in  Figure  5. 
Panels  a,  c,  and  e  contain  the  scaling  functions;  the  wavelets  are  shown  in  the  other  panels. 
Panels  a  and  b  are  the  Haar  wavelets,  and  the  other  panels  contain  the  Daubechies  wavelets. 
They  are  denoted  Da  and  Dio,  and  their  values  are  tabulated  in  Daubechies’  seminal  paper 
(Reference  24).  The  subscripts  denote  the  number  of  vanishing  moments;  two  times  it  is 
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the  minimal  number  of  coefficients  that  are  necessary  to  implement  the  respective 
Daubechies  wavelets  via  the  discrete  wavelet  transform.  Incidentally,  D\  is  the  Haar 
wavelet. 

Figures  6,  7,  8,  and  9  illustrate  the  different  values  of  the  wavelet  transform  operating 
on  the  Bu(i)  realizations.  The  parameter  d  in  the  ARIMA(0,d,0)  model  equals  0.4  in 

Figures  6  and  8  and  -0.4  in  Figures  7  and  9.  Panels  b,  d,  and  f  in  each  of  these  figures  are 
plots  of  the  detail  coefficients  d,(n)  j  =  1,2,3  for  the  corresponding  wavelet  transform.  In 

each  figure  panel  c  is  the  correlation  evaluated  using  the  third  scale  detail  coefficients 
dAn).  Panel  c  in  each  of  these  figures  indicates  that  the  wavelet  transform  does  very'  well 
in  decorrelating  or  whitening  the  fBm  process,  but  it  is  clearly  dependent  on  the  sign  of  the 
parameter  d.  The  ability  of  wavelets  to  favorably  transform  the  fBm  is  the  basis  for  the 
filtering  theory  presented  by  Chou  (Reference  38).  The  reason  for  this  fortuitous 
transformation  is  given  by  Fiandrin  (Reference  23)  and  Tewfik  and  Kim  (Reference  31). 

The  decay  of  the  correlation  coefficients  is  directly  controlled  by  the  vanishing  moment 
properties  of  the  wavelet; 

]tmys(t)dt  =  0  m  =  0,1,-  -,M  - 1  forM  >0 

Panel  e  in  each  figure  is  a  plot  of  the  approximation  coefficients  a;(n)  at  the  scale  j  =  3. 
The  number  of  vanishing  moments  for  D«  and  Dio  are  4  and  10. 

As  a  proposed  model  for  the  detector  noise  on  staring  arrays.  fBm  has  three  parameters 
(the  mean  fi,  the  variance  a ,  and  the  Hurst  coefficient  H)  that  must  be  estimated  directly 
from  measurements.  Of  the  three  parameters,  the  index  H  is  the  most  difficult  to  estimate 
and  is  clearly  the  most  important.  In  Figures  10  and  1 1  the  index  or  Hurst  coefficient  H  is 
estimated  using  the  rescaled  range  and  the  wavelet  transform.  In  panels  a,  c,  and  e,  H  is 
estimated  from  the  slope  on  the  log-log  plot  In  all  panels  a,  c,  and  e,  the  solid  line  connects 
the  values  of  the  variance  at  the  respective  scales  ranging  from  1  to  6  or  7,  and  the  dotted 
lines  are  the  least  squares  fit.  The  vertical  axis  is  the  variance  computed  from  wavelet  detail 
coefficients  and  the  horizontal  axis  is  the  scale.  The  estimate  depends  on  the  equation 
presented  by  Mallat  (Reference  27). 

log2(var(dy(n))  =  (2//)j  +  constant 

This  procedure  is  derived  by  Fiandrin  (Reference  22)  and  is  based  on  the  power  law 
behavior  for  the  wavelet  coefficients  of  fBm,  namely. 
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v*r(dj(n))  =  (°2/{)V(p(HX2j)2H  +  l 

where  VP(H)  is  a  constant  that  depends  on  the  wavelet  and  the  index  H. 

The  simulated  value  of  d  is  0.4,  which  yields  H  =  0.8  based  on  the 
equation  H  =  d  +  1/2. 


In  Figures  10  and  11,  panels  b,  d,  and  f  are  the  log-log  plot  containing  the  rescaled 
range  statistic  %  versus  the  log  of  the  block  size  parameter  x.  H  is  obtained  by 

/  J 


estimating  the  slope.  In  the  previous  section  the  justification  for  the  use  of  the  rescaled 

range  R/S  as  a  statistic  is  given  (i.e.,  %  m  |fj  )•  Although  a  careful  asymptotic  analysis  of 

/  0 
A 

these  two  statistics  for  H  is  apparently  unavailable  in  the  literature,  the  variability  of  these 
estimates  as  a  function  of  sample  size  is  clearly  evident.  Based  on  these  studies  a  sample 
size  of  3,000  or  more  points  is  clearly  desirable  to  help  stabilize  the  estimates.  The  sample 
size  requirement  is  not  an  inordinate  burden  for  infrared  cameras  that  sample  at  a  30  (60) 
Hertz  rate,  because  3,000  samples  can  be  obtained  in  100  (50)  seconds — a  reasonable 
amount  of  time  for  many  applications!  However,  more  efficient  statistical  estimates  with 
confidence  intervals  are  clearly  needed. 


A 

The  wavelet  estimates  for  the  index  H  in  Figures  10  and  1 1  should  be  compared  with 
the  wavelet  estimates  in  Figure  12.  In  all  panels,  the  solid  line  connects  the  values  of  the 
variance  at  the  respective  scales  ranging  from  1  to  7,  and  the  dotted  lines  are  the  least 
squares  fit.  Panels  c,  d,  e,  and  f  use  the  single  wavelet  Dio  at  all  scales.  The  inflection 
point  located  at  scale  2  in  panels  c,  d,  e,  and  f  is  also  observed  in  the  plots  presented  by 
Flandrin  (Reference  22).  He  also  discusses  the  rationale  for  this  bend;  the  interested  reader 
is  referred  to  his  paper.  Panels  b,  d,  and  f  use  each  member  of  the  Daubechies  wavelet 
family,  £>io.  Da,  and  D2  in  sequential  order  for  two  scales.  The  reduction  in  variability 
over  the  Dio  only  transforms  is  clearly  evident.  Again,  the  rational  for  the  reduction  in 
variability  is  discussed  by  Tewfik  and  Kim  (Reference  31).  As  the  scales  progress  from 
fine  to  coarse  the  support  of  the  wavelets  needs  to  be  adjusted  to  match  the  decrease  in 
sample  size  at  each  scale. 

Because  any  measurement  obtained  from  an  infrared  camera  detector  will  include  both 
the  1/f-type  noise  depending  on  the  detector  material  properties  and  other  measurement- 
induced  noise  (e.g.,  analog-to-digital  conversion,  quantizing  effects),  it  is  reasonable  to 
expect  the  periodogram  to  be  corrupted  with  these  noises.  If  the  a  priori  knowledge  of  the 
spectral  density  of  the  measurement  noise  is  unknown,  then  the  problem  becomes  one  of 
jointly  measuring  signal  and  noise.  Womell  does  address  the  problem  by  postulating  an 
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additive  model  of  a  1/f  signal  corrupted  by  stationary  white  Gaussian  noise  (Reference  12). 
In  Figure  13,  the  simulated  additive  model  sample  paths  are  illustrated  in  panels  a  and  c 
with  corresponding  d  values  of  0.4  and  -0.4.  The  noise  variance  is  equal  to  4  in  panel  a  and 
0.0036  in  panel  c.  Panels  a,  b,  c,  and  d  in  Figure  12  should  be  compared  with  the  panels  a, 
b,  e,  and  f  in  Figure  2.  The  main  effect  of  the  noise  is  to  alter  the  rolloff  rate  in  the  high- 
frequency  interval  in  the  periodogram.  The  accumulated  variance  obtained  from  the  sample 
paths  BfjO)  in  panels  a  and  c  is  plotted  in  panels  e  and  f.  The  increase  in  the  variance  is 

another  characteristic  of  persistent  and  antipersistent  behavior. 

Besides  the  index  H,  the  other  two  parameters  in  the  fractional  Brownian  model  are  the 
mean  fJ.  and  the  variance  a.  These  parameters  will  completely  characterize  fBm  provided 
the  marginal  distributions  of  the  increments  (r)  are  Gaussian. 

While  the  marginal  distributions  of  the  measured  Amber  camera  data  will  probably 
never  exactly  match  this  standard,  any  respectable  model  ide  ntification  procedure  should 
define  the  degree  of  approximation. 

Fortunately  some  graphical  techniques  exist  to  examine  the  marginal  distributions. 
Figures  14  and  15  exhibit  three  methods  found  in  S-PLUS.  In  Figure  14  the  first  four 
panels  contain  Q-Q  plots.  A  Q-Q  plot  can  be  used  to  compare  the  degree  of  agreement 
between  two  empirical  distributions  or  can  be  used  to  compare  the  empirical  quantiles  with 
the  quantiles  from  an  ideal  distribution.  A  Q-Q  plot  is  a  plot  of  the  ordered  data  y,  from  the 
sample  {*<>,••  sxat-i}  versus  the  normal  quantiles  yp =  where 

P<  =  0  -^)/n,  i  =  l,2,*",N,  and  O'1  is  the  inverse  of  the  standard  normal  distribution 
(Reference  37).  If  the  shape  of  the  marginal  distributions  for  the  increments  VB^(r)  is 

approximately  normal,  even  in  the  tails,  then  the  empirical  quantile  sample  values  will 
approximate  the  normal  line.  The  Q-Q  plots  for  the  increments  in  panels  a  and  c  are  linear 
even  in  the  tails,  so  the  Gaussian  hypothesis  for  the  marginal  distribution  as  simulated  by 
the  Hosking  ARIMA(0,d,0)  model  is  reasonable.  Clearly,  the  Q-Q  plots  for  the  sample 
paths  are  quite  different  with  the  persistent  model  exhibiting  a  radical  departure  from 
normality. 

Panels  e  and  f  illustrate  how  the  variance  influences  the  range  or  excursion  of  the 
sample  path: 

R{ T)=  max3((f,T)  -  min  X(t, x) 
l<tsr  IStST 
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The  Hosking  fGn  ARIMA(0,d,0)  model  does  not  constrain  the  variance  to  any  particular 
value  in  the  Gaussian  generator,  and  for  modeling  purposes  it  is  useful  to  have  the  variance 
as  an  extra  degree  of  freedom. 

Figure  15  shows  the  histograms  and  density  estimates  for  the  increments  and  sample 
paths.  The  agreement  among  the  three  estimates  of  distribution  shape  is  quite  good.  Of 
course,  any  of  these  exploratory  tools  that  are  used  to  measure  the  shape  are  not  a  substitute 
for  other  tests  for  normality. 

Figures  16  and  17,  respectively,  exhibit  correlation  functions  between  two  different 
realizations  of  fBm  and  their  associated  increments.  In  Figure  16  panels  c  and  f  are  the 
autocorrelations  for  the  corresponding  fBm  sample  paths  in  panels  a  and  b.  The  positive 
coefficients  and  slow  decay  of  the  autocorrelation  functions  is  one  of  the  distinctive 
characteristics  of  processes  with  a  spectral  density  like  ARIMA(O.d.O)  d  <  0.  Most 
traditional  time-series  analysis  has  been  concerned  with  the  property  that  observations 
separated  by  a  long  time  span  are  nearly  independent.  These  types  of  processes  have 
correlation  functions  that  are  typically  like  those  in  Figure  17,  panels  c  and  d.  Panels  d  and 
e  are  the  corresponding  cross  correlations  (with  negative  and  positive  lags)  between  the 
fBm  in  panels  a  and  b.  Panels  d  and  e  basically  illustrate  a  small  but  persistent  degree  of 
long-term  negative  correlation,  which  is  consistent  with  the  apparent  degree  of  cyclical 
behavior  in  the  two  sample  path'.  By  contrast  the  cross-correlation  plots  for  the  increments 
indicate  a  total  lack  of  correlation  between  the  increments. 

Figure  18  illustrates  the  analysis  and  synthesis  processing  versatility  of  the  wavelet 
transform  (see  Figure  1).  Panel  a  is  the  original  signal  fBm  sample  derived  from  a  Hosking 
ARIMA(0,-0.4,0)  model;  panels  b,  d,  and  f  are  the  detail  coefficients  d,(n)  j  =  1,2,3;  and 
panels  c  and  e  are  the  approximation  coefficients  dj(n)  j  =  2,3.  Panels  g  and  h  represent 
the  resulting  signal  after  wavelet  synthesis,  except  that  only  selected  scales  from  either  the 
approximation  coefficients  aj(n)  or  the  detail  coefficients  d;(n)  are  used.  Panel  g  is  the 
signal  that  results  from  reconstructing  from  scale  e;  i.e.,  only  the  coefficients  d,(n)  j  =  3 
are  used  for  synthesis.  Note  how  the  density  of  noise  in  panel  g  is  reduced  from  that  in 
panel  a.  Panel  h  is  the  signal  that  results  from  reconstructing  from  scales  b,  d,  and  f;  i.e., 
only  the  coefficients  dj(n)  j  =  1,2,3,  are  used  for  synthesis.  As  expected,  now  the 
reconstructed  signal  h  only  has  noise — the  low-frequency  drift  has  been  filtered  out  of  the 
signal. 
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AMBER  CAMERA  TEST  RESULTS 


In  this  section  some  pixel  data  taken  from  an  Amber3  AE  4128  imaging  system  are 
tested  to  determine  whether  the  fBm  hypothesis  is  reasonable.  The  Amber  camera  is  a  128 
by  128  indium  antiomide  (InSb)  staring  array  that  operates  at  midwave  infrared 
wavelengths  of  3  to  5.5  /im.  The  imaging  system  includes  an  FPA,  a  transimpedance 
amplifier  (TLA),  a  pour-filled  dewar,  an  f3  lens,  and  an  electronics  unit  for  camera  control. 

Before  the  test  the  vacuum  was  pumped  down  to  10~5  torr.  The  camera  lens  was 
located  about  1  inch  from  the  blackbody,  which  was  set  at  a  constant  temperature  of  20°C 
during  the  measurements.  The  blackbody  has  a  4-  by  4-inch  surface  area  and  is  inv  rfaced 
to  an  EO  Industries  blackbody  temperature  controller.  The  blackbody  and  temperature 
controller  are  accurate  to  ±0.01  °C.  Two  sets  of  measurements  were  conducted.  Each 
distinct  set  was  quantized  to  12  bits  of  sampled  data  from  the  array.  The  first  measurements 
consisted  of  two  different  sets  of  image  sequences,  namely,  a  burst  set  that  contained 
1,024  images  with  a  sampling  period  of  1.4649  seconds,  and  long-term  measurement  set 
that  contained  800  images  with  a  sampling  period  of  60.6  seconds.  The  burst  data  was 
taken  prior  to  each  extended  measurement  The  dewar  was  filled  with  liquid  nitrogen  at  the 
beginning  of  each  test  and  was  not  refilled  during  the  13.5  hours  of  each  test  This 
extended  period  without  any  refill  is  within  the  Amber  guidelines  for  camera  operation; 
however,  saturation  was  clearly  evident  well  beyond  the  800  images  retained  in  the  long¬ 
term  measurements  database. 

The  second  set  of  measurements  required  modifying  the  digital  memory  buffers  so  that 
the  transfer  rate  from  the  array  elements  to  the  storage  memory  was  commensurate  with  the 
maximum  sampling  rate  of  the  array,  which  is  54  Hertz.  Because  the  original  memory 
buffers  that  were  used  for  the  data  transfer  from  the  array  were  designed  for  another  task, 
to  achieve  the  54  Hertz  rate  a  compromise  solution  was  designed  by  Douglas  Lamb  of  the 
Advanced  Signal  Processing  Branch  at  NAWCWPNS  China  Lake.  He  was  able  to  sample 
a  16-  by  16-element  subairay  at  the  maximum  rate.  This  data  set,  called  the  54  Hertz  set,  is 
discussed  in  the  final  paragraphs  of  this  section. 

Table  3  lists  the  key  camera  parameter  settings  that  were  used  during  all  of  the  tests. 

In  both  Figures  19  and  20,  the  panels  labeled  “camera  data”  are  the  time-sampled  output 
from  an  individual  pixel  on  the  Amber  camera  array.  In  each  of  these  plots,  the  horizontal 
axis  represents  the  number  of  samples,  and  the  vertical  axis  is  a  gray  scale. 


3  AMBER  Engineering,  Inc.,  5756  Thomwood  Drive.  Goleta,  California,  93 1 17-3802. 
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TABLE  3.  Camera  Parameter  Settings. 


Camera  parameters 

Setting 

T1A  gain 

3 

TLA  offset 

0 

Integration  time 

32 

Frame  rate 

54 

Global  gain 

1 

Global  offset 

0 

The  purpose  of  the  panel  plots  in  Figure  19  is  to  illustrate  that  the  estimate  of  the  Hurst 
index  H  can  be  extremely  sensitive  to  the  number  of  samples.  Panel  a  is  a  fBm  with  a 
Hosking  ARIMA(0,0.3,0)  model,  and  panel  b  is  the  periodogram  with  the  dotted  line 
representing  least  squares  fit  for  the  estimate  of  H.  The  estimate  is  in  good  agreement  with 
the  proper  value.  Panel  c  is  a  segment  of  the  fBm  in  panel  a  that  is  selected  for  its  bow, 
which  is  a  good  approximation  to  the  burst  data.  Panels  d  and  f  are  the  associated 
periodograms,  both  of  which  give  a  value  for  the  index  H  that  differs  from  the  value  in 
panel  b.  Of  course,  the  estimate  for  H  in  panel  d  is  clearly  wrong!  By  way  of  comparison 
the  fraediff  program  with  the  panel  c  sample  was  able  to  estimate  H  correctly.  Three 
lessons  are  illustrated  in  Figure  19.  First,  sample  size  is  important.  Second,  the 
periodogram  may  not  always  be  a  reliable  estimator  of  the  index  H,  especially  if  the 
sampled  data  set  is  too  small.  Third,  the  noise  does  influence  the  high-frequency  values  of 
the  periodogram  in  a  manner  similar  to  that  predicted  in  Figure  13.  Unfortunately,  no  clear 
guidelines  exist  about  these  matters,  other  than  basic  Fourier  estimates  for  the  maximum 
and  minimum  frequency  band  that  must  be  measured  to  obtain  the  best  estimate  of  the  index 
H .  Obviously,  the  lowest  frequency /min  is  determined  by  1/T  where  T  is  the  total  sample 
time,  so  a  larger  value  of  T  means  a  lower  frequency.  The  first  two  entries  in  Table  4  are 
the  parameters  for  the  sampled  data  sets,  and  the  last  two  entries  are  the  minimum  and 
maximum  frequencies  in  the  measured  frequency  band.  Here  A/  is  the  sampling  interval. 
The  first  two  rows  represent  the  burst  and  long-term  data  sets,  respectively,  while  the  last 
three  rows  are  the  parametric  values.  A  60-Hertz  sampling  rate  is  assumed. 

TABLE  4.  Fourier  Sampling  Parameters. 


Number  of  samples 

Ar,  seconds 

/min  •  Hertz 

/max  •  Hertz 

1024 

1.4649 

0.00066664 

0.68264 

800 

60.6 

0.00002062706 

0.0165 

3,600 

0.016667 

0.016667 

60 

7,200 

0.016667 

0.008333 

60 

360,0000 

0.016667 

0.000167 

60 
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In  Figure  20,  the  variability  of  the  statistic  H  due  to  data  smoothing  is  illustrated.  These 
changes  can  be  anticipated  by  consulting  Figure  2. 

In  Figure  21  panels  a,  c,  and  e  are  plots  of  640  time-sampled  values  from  pixels  (10,9), 

(10,10),  and  (10,11),  respectively,  taken  from  the  long-term  measurements.  Panels  b,  d, 
and  f  are  the  periodograms,  and  the  index  H  is  the  slope  of  the  dotted  line  obtained  by  a 
least  squares  fit  The  agreement  among  the  indices  is  consistent  with  the  hypothesis  that  the 
spatial  noise  on  the  Amber  camera  is  fBm  with  a  mean  power  law  of  Table  5 

includes  estimates  of  the  index  H  taken  over  a  larger  sample  of  pixels. 

In  Figure  22  panels  a,  c,  and  e  are  plots  of  640  time-sampled  values  from  pixels  (10,9), 

(10.10) ,  and  (10,11),  respectively,  taken  from  the  long-term  measurements. 

Panels  b,  d,  and  f  are  the  Q-Q  plots  for  the  increments.  The  spacing  of  the  bars  on  in 
the  Q-Q  plots  is  caused  partially  by  the  coarseness  of  the  quantization  for  the  measured 
interval  (i.e.,  the  12-bit  quantization  is  spread  over  a  larger  range).  Clearly  the  increments 
are  well  approximated  by  the  Gaussian  shape.  The  corresponding  histograms  and  density 
estimates  are  shown  in  Figure  23  and  clearly  confirm  the  quantile  plots  prediction. 

Panels  a  and  c  in  Figure  24  are  plots  of  640  time-sampled  values  from  pixel  (10,10) 
taken  from  the  long-term  measurements.  The  other  panels  show  the  reasonable  agreement 
among  the  Hurst,  wavelet,  and  periodogram  estimates  for  the  index  H .  The  new  addition 
is  panel  f,  which  is  a  plot  of  the  correlation  function  for  the  increments.  The  Hosking  fGn 
ARIMA(0,0.d,0)  model  predicts  that  the  correlation  coefficients  will  all  be  positive  for  d  > 
0  and  for  d  <  0,  and  the  coefficient  will  be  negative  and  the  follow  the  positive  axis  as  a 
function  of  the  lag.  The  Amber  camera  correlation  is  consistent  with  the  antipersistent  case, 
namely  (d  <  0)  and  should  be  compared  with  Figure  17,  panels  e  and  f.  The  agreement 
between  the  simulated  fBm  and  the  Amber  camera  in  this  case  is  excellent. 

Figures  25  and  26  exhibit,  respectively,  the  correlation  functions  between  pixels 

(10.10)  and  (10,11)  and  their  associated  increments.  In  Figure  25  panels  a  and  f  are  the 
autocorrelations  for  the  spatial  noise  in  panels  c  and  b.  The  positive  coefficients  and  slow 
decay  of  the  autocorrelation  functions  is  one  of  the  distinctive  characteristics  of  processes 
that  the  fBm  simulation  clearly  exhibited  in  Figures  3  and  16.  Thus,  the  observations 
separated  by  a  long  time  span  are  more  persistent  and  are  clearly  not  negligible.  Time  series 
with  the  slow  decay  exhibit  cycles  and  changes  of  level  of  all  orders  of  magnitude — a 
description  that  matches  both  the  Amber  camera  data  and  the  fBm.  Panels  d  and  e,  the 
corresponding  cross  correlations  between  panels  a  and  b,  illustrate  a  persistent  degree  of 
long-term  positive  correlation,  which  is  consistent  with  the  apparent  degree  of  cyclical 
behavior  in  the  two  sample  paths. 
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TABLE  5.  Amber  Camera  Statistics. 


Pixel  s 

ample  path 

statistics 

1!  Increments  statistics 

_ 1*12, . 

Min 

Max 

Mean 

Var 

|  Min 

Max 

Var 

it 

(»,») 

983 

1017 

998 

44  43 

if  ■ 9 

8 

8.02 

0  3886 

(9.10) 

823 

868 

843 

71.44 

1  -13 

1  4 

24  65 

0.3324 

(9,11) 

942 

97* 

959 

52  45 

1  0 

1  12 

0  37n 

(10.9) 

959 

992 

973 

47.32 

8 

8.06 

0.3594 

(10,10) 

844 

878 

859 

45.99 

.  1  m 

1  1 

9  98 

0  3579 

(10.11) 

806 

*41 

819 

51.11 

■  ■  v„ 

8  31 

0  4035 

(11.9) 

1104 

1141 

1 1  19 

52.12 

M 

■  ' 

8.59 

0  3723 

(11.10) 

856 

890 

869 

47.11 

H  m 

9  94 

0  3798 

948 

9*3 

965 

50.36 

•so 

1  10 

9  62 

0  3612 

889 

922 

904 

46.23 

-  9 

7  94 

0.3779 

733 

773 

750 

52.68 

*13 

1  8 

1?  02 

0  2975 

813 

848 

82? 

44.88 

*9 

9 

9  17 

0  332 

(21,22) 

*17 

*53 

831 

41  55 

-10 

9 

8  96 

0.4241 

(22,21) 

910 

943 

924 

49.33 

*  9 

9 

8  42 

0  4252 

(32.32) 

912 

946 

928 

50.0* 

-1  1 

1  1 

9.52 

0  3994 

(20.64) 

899 

935 

915 

50.15 

9 

8.23 

0  371 

(21.63) 

*77 

910 

893 

45.05 

10  86 

0  3749 

(21.64) 

938 

973 

954 

55  41 

9  34 

0  3531 

(21.65) 

1049 

1085 

1064 

54  3* 

11. 03 

0  35*5 

(22.64) 

850 

890 

868 

51.12 

9.33 

0  3722 

(32.64) 

909 

950 

92* 

56.93 

mm 

1  3 

16.38 

0  3477 

(10.110) 

496 

536 

511 

SI  33 

1 

1  1 

11.17 

0  3897 

(20.105) 

657 

693 

671 

43.01 

*t  i 

9 

10.12 

0.3478 

(21.104) 

577 

6  10 

593 

42.93 

-ii 

1  1 

10.98 

0.3482 

(21.105) 

602 

637 

615 

39.07 

-9 

9 

8.29 

0  3263 

(21,106) 

419 

456 

436 

44.35 

-10 

9 

9.32 

0  3355 

(22.105) 

673 

70* 

6*6 

47.03 

-1 

7 

7.8 

0  3355 

(32.95) 

944 

979 

959 

50.13 

-11 

12 

n.75 

0  3437 

(63.21) 

1301 

1346 

1311 

10.42 

-1 1 

12 

14.52 

0  3472 

(64.20) 

991 

102* 

1007 

55.25 

•13 

1  1 

9.93 

0  3573 

(64,22) 

1025 

1059 

1040 

50.31 

-10 

10 

9.25 

0.3521 

(64.32) 

1105 

1143 

1122 

41.96 

•9 

10 

10.41 

0.3907 

(65.21) 

1213 

1249 

1229 

52.7* 

-* 

13 

*.22 

0  3*27 

(63.64) 

1110 

1150 

1125 

49.16 

-9 

10 

8  4 

0.3766 

(64.64) 

1222 

1256 

1237 

41.49 

-10 

10 

9.82 

0.389 

(64.63) 

12*9 

1330 

1306 

57.(8 

-17 

14 

22.91 

0.3175 

(64,65) 

1316 

1333 

1333 
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Correlation  functions  that  are  typically  like  those  in  Figures  26  panels  c  and  d  mean  that 
observations  separated  by  a  long  time  span  are  nearly  independent.  The  cross  correlation 
plots  for  the  increments  indicate  a  total  lack  of  correlation  between  the  increments.  These 
panels  again  show  a  remarkable  agreement  with  the  correlation  functions  increments  for 
fBm  plotted  in  Figure  17  in  the  same  panels. 

Table  5  contains  the  summary  statistics  from  63  Amber  camera  pixels  from  the  128  x 
128  array.  Each  pixel  sample  is  a  discrete  time  sequence  with  640  samples.  The  First 
column  in  the  table  represents  the  pixel  coordinates  with  indices 

[a,.,]  i  =  0,...,  127, j  =  0,...,127 

The  statistics  are  obtained  for  each  pixel  with  its  discrete  time  sequence  of  640  samples. 
The  statistics  for  the  columns  labeled  “sample  path”  are  for  these  sequences,  and  the 
increments  statistics  are  obtained  after  differencing  the  time  sequences  by  applying  the 
difference  operator  defined  in  the  FBm,  Wavelets,  and  ARIMA  Models  section.  The 
columns  labeled  min,  max,  mean,  and  var  are  the  minimum  values,  the  maximum  values, 
the  sample  mean,  and  the  sample  variance  for  each  discrete  pixel  time  sequence  and  the 
corresponding  increments.  The  column  labeled  H  is  the  fBm  index  H  estimate  discrete  time 
sequence. 

Table  6  presents  the  statistics  for  each  of  the  columns  in  Table  5.  As  is  evident  the 
index  H  has  a  mean  value  of  0.3660,  yielding  a  exponent  of  1.732  for  the  fBm  spectrum. 
The  mean  value  of  the  variance  for  the  fBm  process  is  10.52.  It  is  more  variable  because 
some  of  the  pixels  have  some  outliers,  as  an  examination  of  the  maximum  and  minimum 
values  suggests.  Of  course  the  outliers’  presence  should  be  checked  by  computing  the  Q-Q 
plots.  Recall  that  the  variance  for  fBm  governs  the  marginal  distribution  for  the  increments. 

TABLE  6.  Amber  Camera  Global  ^atistics. 


Pixel  sample  path  statistics  _ ][_ _ Increments  statistics 


Min 

Max 

Mean 

Var 

Min 

Max 

Var 

H 

Mean 

1014.5 

1050.8 

1030.38 

51.19 

-10.1 

10.13 

10.52 

0.366 

Variance 

46005.6 

46140 

46206.7 

36.98 

3.23 

3.89 

11.11 

0.0009 

StdDev 

214.49 

214.8 

214.96 

6.08 

1.8 

1.97 

3.33 

0.0293 

Min 

419 

456 

436 

39.07 

-17 

7 

7.39 

0.2975 

Max 

1351 

1387 

1367 

80.42 

-8 

18 

24.65 

0.4362 

In  Figure  27  the  panels  a,  c,  and  e  are  the  time-sampled  output  from  an  individual  pixel 
on  the  Amber  camera  array.  In  each  of  these  plots  the  horizontal  axis  represents  the  number 
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of  samples,  and  the  vertical  axis  is  a  gray  scale.  All  of  the  data  in  Figures  27,  28,  29,  30, 
and  31  were  sampled  at  54  Hertz. 

In  Figure  27  panels  a,  c,  and  e  are  time  samples  from  the  Amber  camera  data  from  three 
adjacent  pixels  on  the  same  row.  Panels  b,  d,  and  f  are  the  periodograms,  and  the  index  H 
is  the  slope  of  the  dotted  line  obtained  by  a  least  squares  fit.  The  agreement  among  the 
indices  is  consistent  with  the  hypothesis  that  the  spatial  noise  on  the  Amber  camera  is  fBm. 
The  estimates  for  H  are  0.385304,  0.158828,  and  0.356435,  respectively,  for  pixels 
(8,7),  (8,8)  and  (8,9).  Using  the  formula  a  =  2  H  +  1  yields  a  power  spectrum  inversely 
proportional  to  the  frequency  f~a  with  values  1.7706,  1.3176,  and  1.7128,  respectively, 
for  pixels  (8,7),  (8,8)  and  (8,9).  Note  that  pixel  (8,8)  exhibits  some  degree  of  flicker 
noise. 

Panels  b,  d,  and  f  in  Figure  28  are  the  Q-Q  plots  for  the  increments.  The  spacing  of  the 
bars  on  in  the  Q-Q  plots  is  partially  due  to  the  coarseness  of  the  quantization  for  the 
measured  interval  (i.e.,  the  12-bit  quantization  is  spread  over  a  larger  range).  Clearly  the 
increments  arc  well  approximated  by  the  Gaussian  shape.  The  corresponding  histograms 
and  density  estimates,  shown  in  Figure  29,  clearly  confirm  the  quantile  plots  prediction. 

As  a  confirmation  of  the  predicted  behavior  of  wavelet  processing  on  fBm,  Figures  30 
and  31  illustrate  the  wavelet  transform  ability  to  “  decorrelate  or  whiten”  fBm  as  predicted 
by  Flandrin  (Reference  22).  Panels  b,  d,  and  f  in  each  of  these  figures  are  plots  of  the  detail 
coefficients  dj(n)  j  =  1,2,3  for  the  corresponding  wavelet  transform.  In  each  figure  panel  c 
is  the  correlation  evaluated  using  the  third  scale  detail  coefficients  di(n).  Panel  c  in  each  of 
these  figures  compares  very  favorably  with  the  simulated  counterparts  in  Figures  6  and  8. 
These  results  clearly  support  the  Flandrin ’s  claim. 


CONCLUSIONS 


The  most  important  contribution  of  this  report  is  the  demonstration  that  Daubechies 
wavelets  do  filter  the  fixed  pattern  noise  in  a  predictable  fashion.  This  claim  is  based  on  the 
results  obtained  by  applying  her  wavelet  transform  to  the  temporal  fixed  pattern  noise 
measurements  taken  by  an  InSb  Amber  FPA.  The  1/f  character  of  fixed  pattern  temporal 
noise  measured  in  the  3-  to  5.5- fim  region  by  a  128  by  128  InSb  staring  array  has  been 
identified  as  fractional  Brownian  motion.  Wavelet  processing  successfully  decomelated  the 
fixed  pattern  noise.  The  results  of  this  report  can  also  serve  as  basis  for  fixed  pattern  noise 
models  that  can  be  used  in  missile  simulations  to  test  various  detection  and  nonuniformity 
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compensation  algorithms.  All  three  sets  of  data  were  measured  during  this  study,  and  each 
exhibits  the  presence  of  long-term  correlation,  which  is  a  necessary  condition,  but  not 
sufficient  for  having  a  1/f-type  spectrum.  The  correlation  of  the  three  data  sets  is  also 
predictably  changed  after  the  sets  are  wavelet-transformed.  Again,  the  change  in  correlation 
is  a  necessary  condition  and  is  another  confirmation  of  the  fBm  model. 

An  important  implication  of  the  wavelet  characterization  of  fixed  pattern  noise  is  that  it 
suggests  an  optimum  detection  strategy  in  imaging  arrays.  Presently,  the  optimum  detection 
strategy  in  the  presence  of  fixed  pattern  noise  and  clutter  is  a  key  unresolved  problem. 
Although  conventional  detection  algorithms  will  suffice  for  the  white  noise  and  shot  noise 
present  on  these  arrays,  the  fixed  pat^.m  noise  requires  novel  techniques  for  optimum 
detection.  The  rationale  for  the  detection  >f  a  point  target  in  clutter  proceeds  as  follows  and 
is  suggested  by  the  radar  methodology  tor  this  problem,  which  is  well  known  and  which, 
by  analogy,  suggests  the  following  steps  for  infrared  arrays.  First,  solve  the  detection 
problem  when  only  fixed  pattern  noise  is  present  by  characterizing  the  noise  on  the  array. 
Next,  define  an  optimum  detection  algorithm.  This  paradigm  is  exactly  parallel  to  standard 
processing  that  has  been  successfully  used  in  radar-detector  design.  Next,  formulate  an 
adaptive  detection  algorithm  recognizing  the  array  noise  limits.  Then  maximize  the  detection 
signal-to-clutter  ratio. 

This  report  provides  a  unified  review  of  fractional  Brownian  motion,  wavelets,  and  a 
low-frequency  spectrum  1/f-type  noise  of  the  form  co~a  1  <  a  <  2 .  By  comparison,  most 
traditional  noise  processes  are  characterized  by  short-term  correlation  with  rapid  decay, 
whereas  a  noise  sequence  with  a  1/f-type  spectrum  has  a  correlation  function  that  decays 
very  slowly  and  has  nonzero  values  at  long  distances.  It  is  our  contention  that  the  1/f-type 
noise  is  a  major  contributor  to  array  nonuniformity.  Moreover,  the  long-term  correlation 
structure  of  1/f-type  noise  precludes  the  direct  application  of  conventional  point  target 
detection  algorithms. 

The  report  characterizes  the  1/f-type  noise  as  fractional  Brownian  motion.  Both  the 
theory  and  the  examples  presented  in  the  report  show  that  wavelet  transforms  can  transform 
a  low-frequency,  1/f-type  noise  spectrum  of  the  form  co~a  l<  a  <2  into  a  flat  white- 
noise-type  spectrum.  In  other  words,  the  1/f-type  noise  has  been  decorrelated. 
Decorrelation  is  important  because  it  means  that  the  wavelet  transform  acts  as  a  simple  1/f 
filter.  The  demonstration  that  the  wavelet  transform  decorrelates  the  fixed  pattern  noise 
means  that  the  array  nonuniformity  has  been  reduced  and  standard  detection  algorithms  for 
decorrelated  white  noise  can  be  applied  after  the  wavelet  transform.  The  implications  for 
point  target  detection  in  fixed  pattern  noise  means  that  the  noise  on  all  conventional  staring 
arrays,  whether  they  be  InSb,  HgCdTe,  or  PtSi,  can  be  decorrelated.  Thus,  the  apparent 
lack  of  1/f-type  noise  in  the  PtSi  arrays  would  not  be  such  an  obvious  major  advantage, 
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because  the  1/f-type  noise  on  InSb  and  HgCdTe  can  be  converted  to  white  noise  by  wavelet 
processing. 

The  fractional  Brownian  motion  model  links  the  covariance  of  the  sample  path  or  time 
history  as  it  evolves  in  time  with  its  power  spectrum.  It  is  a  parasimonious  model  (it 
depends  on  two  parameters)  for  1/f  noise.  In  fact,  one  important  feature  of  fBm  is  that  it 
provides  a  model  of  how  long-range  correlation  is  directly  related  to  the  parameters  of  the 
spectrum.  By  comparison,  most  traditional  noise  processes  are  characterized  by  short-term 
correlation  with  rapid  decay,  whereas  a  noise  sequence  with  a  1/f-type  spectrum  has  a 
correlation  function  that  decays  very  slowly  and  has  nonzero  values  for  long  lags.  While 
noise  can  be  classified  according  to  its  power  spectrum,  it  is  important  to  consider  the 
relation  between  the  time  domain  and  the  frequency  domain.  An  important  implication  of 
the  fBm  model  is  that  the  power  spectrum  (representing  the  frequency  domain)  is  linked  via 
the  Fourier  transform  with  the  covariance  (representing  the  time  domain).  This  linkage 
predicts  a  time-dependent  envelope  for  the  maximum  excursion  of  the  sample  path  of  each 
pixel  as  measured  by  the  variance.  This  variance  formula  provides  a  link  between  spatial 
noise  on  the  array  and  the  time- varying  noise  associated  with  an  individual  pixel,  because 
the  model  predicts  that  the  measured  responses  of  the  individual  pixels,  which  are  observed 
by  viewing  a  constant  source,  will  eventually  drift  randomly  away  from  their  original 
nominal  value.  In  fact,  the  power  spectrum  as  measured  by  the  periodogram  of  the 
individual  pixels  provides  an  excellent  indicator  of  the  type  of  spatial  noise  that  will  be 
present  on  the  array,  especially  if  the  spatial  noise  is  measured  by  the  ratio  of  the  sampled 
variance  across  the  array  divided  by  the  sampled  mean  across  the  array.  For  example,  if  the 
spectrum  is  white,  then  this  ratio  for  the  spatial  noise  will  be  statistically  constant,  but  if  the 
spectrum  is  of  the  form  0)~a  with  1  <  a  <  2,  then  the  ratio  will  increase  at  a  rate  that  is 
directly  related  to  the  exponent  a.  In  this  sense  temporal  noise  and  spatial  noise  are  not 
independent! 
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FIGURE  2.  FBms  and  Their  Periodograms. 
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FIGURE  3.  FBms  and  Their  Autocorrelation  Functions. 
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FIGURE  4.  Increments  of  FBms  and  Their  Correlation  Functions. 
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FIGURE  5.  Scale  and  Wavelet  Functions  of  Several  Wavelets. 
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FIGURE  6.  FBm  With  d  =  0.4  and  Haar  Wavelet 
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FIGURE  7.  FBm  With  d  =  -0.4  and  Haar  Wavelet. 
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FIGURE  8.  FBm  With  cl  =  0.4  and  Daubcchies  10  Wavelet. 
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FIGURE  9.  FBm  With  d  =  -0.4  and  Daubechies  10  Wavelet. 


tagqvwtano*)  log^nnanc*)  tog2(vtitonca) 


NAWCWPNSTP8103 


(c)  Wavelet  estimated  H  of  fBm  of  2000  samples  is  0.846639 
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(d)  Hurst  estimated  H  of  fBm  of  2000  samples  is  0.839091 
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(f)  Hurst  estimated  H  of  fBm  of  3000  samples  is  0.823138 


FIGURE  10.  Estimation  of  H  of  fBm  With  d  =  0.4  (H  =  0.9)  by  Hurst  Rescaled 
Range  and  Daubechies  10  Wavelet 
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FIGURE  1 1 .  Estimation  of  H  of  fBm  With  d  =  0.4  (H  =  0.9)  by  Hurst  Rescaled 
Range  and  Daubechies  10  Wavelet. 
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(b)  Mixed  Daubechies  wavelets  asfimatad  H  of  !Bm  of  3000  samples 
is  0.097576 


(d)  Mixed  Daubechtas  wavelets  estimated  H  of  fBm  of  8000  samples 
is  0.108194 


(f)  Mixed  Daubecbies  wavelets  ass  naiad  H  of  fBm  of  8000  samples 
is  0.094145 


FIGURE  12.  Estimation  of  H  of  fBm  With  d  =  -0.4  (H  =  0.1)  by  Daubechies  10  and 
Mixed  Daubechies  Wavelets. 
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FIGURE  13.  Sensitivity  of  fBm  to  Noise. 
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FIGURE  14.  Normal  Q-Q  Plots  and  Variance  Step  Size  Effect. 
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FIGURE  15.  Histograms  and  Estimated  Probability  Density  Functions  of  fBms 
and  Their  Increments. 
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(a)  On*  IBm  with  d  » -0.3  (b)  Another  fBm  with  d  -  -03 
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(e)  Cross  correlation  function  between  (b)  and  (a) 
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FIGURE  16.  Correlation  Functions  Between  fBms  With  d  =  -0.3. 
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FIGURE  17.  Correlation  Functions  Between  Increments  of  fBm  With  d  =  -0.3. 
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(9)  Reconstruction  from  •  scale  (•)  (h)  Reconstruction  from  d  scales  (b),  (d).  and  (0 

FIGURE  18.  Decomposition  and  Reconstruction  of  fBm  Using  Daubechies  10  Wavelet. 
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FIGURE  19.  Comparison  of  the  Estimated  H  for  a  Short-Term  and  a  Long-Term  fBm. 
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I09'0(w) 

(f)  Pariodogram  ol  (a).  Estimaiad  H  Is  0.500936 


FIGURE  20.  Power  Spectntms  for  Amber  Camera  Data  and  Their  Spline  Fits. 
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woe*) 

(a)  Pixal  at  (10,11)  (f)  Pariodogtam  of  (a).  Etfmatad  H  is  0.403508 

FIGURE  21.  Long-Term  Amber  Camera  Data  and  Their  Periodograms. 
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FIGURE  22.  Amber  Camera  Data  and  Their  Normal  Q-Q  Plots. 
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(a)  Hologram  of  increment  of  camera  data,  ptxal  (10,9) 


-to  -5  o  s  10 


(b)  Estimated  probabiftity  density  function  of 
increment  of  camera  data,  pixel  (10.9) 


(c)  Histogram  of  increment  of  camera  data,  pixel  (10,10) 


(d)  Estimated  protnbiltity  density  function  of 
increment  of  camera  data,  pixel  (10.10) 


(f)  Estimated  probabiftity  density  function  of 
Increment  of  camera  debt,  pixel  (10.11) 


FIGURE  23.  Histograms  and  Estimated  Probability  Density  Functions  of  Long-Term 
Amber  Camera  Data. 
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(d)  Variances  of  wavelet  d  scale  of  (a).  Estimated  H  is  0.37606 1 
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(e)  Hurst  coefficient  of  <c).  Estimated  His  0  412115 


FIGURE  24.  Comparison  of  the  Estimated  H  for  Long-Term  Amber  Camera  Data. 
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(•)  Crets  corrafaton  function  tawwan  (b)  and  (a) 


(f)  Auto  correlation  (unction  of  (b) 


FIGURE  25.  Correlation  Functions  Between  Long-Term  Amber  Camera  Data 
Pixel  (10,10)  and  Pixel  (10,11). 
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FIGURE  26.  Correlation  Functions  Between  Long-Term  Amber  Camera  Data 
Pixel  (10,10)  and  Pixel  (10,1 1). 
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(c)  AMBER  camera  data,  pixel  (B,8) 


logtOiw) 

(d)  Periodogram  ot  (c).  Estimated  H  is  0.158828 
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<a)  AMBER  camera  data,  pixel  (8,7) 
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0  2000  4000  6000 

(c)  AMBER  camera  data,  pixel  (8,8) 


'3  2000  4000  6000 

(e)  AMBER  camera  data,  pixel  (8,9) 


FIGURE  28.  54-Hert7  Amber  Camera  Data  and  Their  Normal  Q-Q  Plots. 
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t 


(b)  Estimated  probability  density  function  ol 
increment  of  camara  data,  pixel  (8.7) 


(d)  Estimated  probabStity  density  lunction  of 
moment  of  camera  data,  pixel  (B.8) 


(f)  Estimated  probability  density  lunction  of 
increment  of  camera  data,  pixel  (8,9) 


FIGURE  29.  Histograms  and  Estimated  Probability  Density  Functions  of 
54-Hertz  Amber  Camera  Data. 
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(8)  Haar  wavelet  a  scale  3  of  (a)  (f)  Haar  wavelet  d  scale  3  ot  (a) 


FIGURE  30.  54-Hertz  Amber  Camera  Data  and  Haar  Wavelet. 
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(f)  Daubechies  10  wavelet  d  scale  3  ol  (a) 


FIGURE  31. 54-Hertz  Amber  Camera  Data  and  Daubechies  10  Wavelet 
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